rule mark_duplicates:
    input:
        sorted_bam="{sample}.sorted.bam",
    output:
        sorted_markdup_bam="{sample}.sorted.markdup.bam",
        sorted_markdup_bai="{sample}.sorted.markdup.bai",
        markdup_metrics="{sample}.markdup.metrics",
    params:
        java_options=config.get("java_options", "")
    shell:
        """
        gatk {params.java_options} MarkDuplicates \
            -I {input.sorted_bam} \
            -O {output.sorted_markdup_bam} \
            -M {output.markdup_metrics} \
            --CREATE_INDEX true
        """

rule base_recalibrator:
    input:
        sorted_markdup_bam="{sample}.sorted.markdup.bam",
        ref_fasta=config["ref_fasta"],
        dbsnp_vcf=config["dbsnp"],
        hapmap_snps_vcf=config["hapmap_snps"],
        omni_snps_vcf=config["omni_snps"],
        onekg_snps_vcf=config["onekg_snps"],
        mills_indels_vcf=config["mills_indels"],
        known_indels_vcf=config["known_indels"],
    output:
        recal_tbl="{sample}.recal.tbl",
    params:
        java_options=config.get("java_options", "")
    shell:
        """
        gatk {params.java_options} BaseRecalibrator \
            -I {input.sorted_markdup_bam} \
            -O {output.recal_tbl} \
            -R {input.ref_fasta} \
            --known-sites {input.dbsnp_vcf} \
            --known-sites {input.hapmap_snps_vcf} \
            --known-sites {input.omni_snps_vcf} \
            --known-sites {input.onekg_snps_vcf} \
            --known-sites {input.mills_indels_vcf} \
            --known-sites {input.known_indels_vcf} \
            -OQ true
        """

rule apply_bqsr:
    input:
        sorted_markdup_bam="{sample}.sorted.markdup.bam",
        ref_fasta=config["ref_fasta"],
        recal_tbl="{sample}.recal.tbl",
    output:
        sorted_markdup_recal_bam="{sample}.sorted.markdup.recal.bam",
        sorted_markdup_recal_bai="{sample}.sorted.markdup.recal.bai",
    params:
        java_options=config.get("java_options", "")
    shell:
        """
        gatk {params.java_options} ApplyBQSR \
            -I {input.sorted_markdup_bam} \
            -O {output.sorted_markdup_recal_bam} \
            -R {input.ref_fasta} \
            --bqsr {input.recal_tbl} \
            --static-quantized-quals 10 \
            --static-quantized-quals 20 \
            --static-quantized-quals 30 \
            --static-quantized-quals 40 \
            --static-quantized-quals 50 \
            -OQ true \
            -OVI true
        """

rule haplotype_caller:
    input:
        sorted_markdup_recal_bam="{sample}.sorted.markdup.recal.bam",
        ref_fasta=config["ref_fasta"],
    params:
        java_options=config.get("java_options", "")
    output:
        gvcf="{sample}.gvcf.gz",
        gvcf_tbi="{sample}.gvcf.gz.tbi",
        call_bam="{sample}.call.bam",
        call_bai="{sample}.call.bai"
    shell:
        """
        gatk {params.java_options} HaplotypeCaller \
            -I {input.sorted_markdup_recal_bam} \
            -O {output.gvcf} \
            -R {input.ref_fasta} \
            --bam-output {output.call_bam} \
            -ERC GVCF \
            -OVI true
        """